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ABSTRACT 

Pulsar timing is a promising technique for detecting low frequency sources of gravitational waves. 
Historically the focus has been on the detection of diffuse stochastic backgrounds, such as those 
formed from the superposition of weak signals from a population of binary black holes. More recently, 
attention has turned to members of the binary population that are nearer and brighter, which stand 
out from the crowd and can be individually resolved. Here we show that the timing data from an array 
of pulsars can be used to recover the physical parameters describing an individual black hole binary 
to good accuracy, even for moderately strong signals. A novel aspect of our analysis is that we include 
the distance to each pulsar as a search parameter, which allows us to utilize the full gravitational 
wave signal. This doubles the signal power, improves the sky location determination by an order of 
magnitude, and allows us to extract the mass and the distance to the black hole binary. 
Subject headings: gravitational waves, pulsars: general 



1. INTRODUCTION 

The detection of gravitational waves via the observa- 
tion of pulse arrival times from pulsars was first proposed 
in the 1970s (lEstabrook fc Wahlqvn^[T975t ISa^hir3tl978l 
iDetweilerl 11979! ) . Recent improvements in pulsar timing 
have made this method one of the best candidates for the 
first detection of gravitational waves. Pulsar timing ar- 
rays (PTAs) are sensitive to gravitational wave frequen- 
cies between / — 10~ 9 Hz to f — 10~ 6 Hz. The lower 
bound is set by the observation time, 30 years being a 
reasonable cut-off, while the upper bound is set by the 
sample rate, with once a week as the goal. There is lit- 
tle motivation to go to higher sample rates as the PTAs 
operate in the short wavelength limit - the distance to 
each pulsar is larger than the wavelength of the gravita- 
tional waves - and the sensitivity falls off as 1/f across 
the band. 

The PTA operating range makes it a excellent 
tool to search for stochastic backgrounds produced 
by a population of slowly evolving supermassive black 
hole binaries ( MBHBs) jR aiagopal & Romani 1995; 
Wvithe fc Loebl 120031: Uaffe fc Backer! 120031: iWen et al l 
20091) . The idea is that the fluctuating time of ar- 
rivals caused by the gravitational waves will produce 
a correlated response across the PTA, with a charac- 
teristic dependence on the angle between each pair of 
pulsars (|Hellings fc Downs! I1983D . Considerable work 
has gone into producing bo unds on the amplitude 
of the stochastic bac kground (iHellings fc Downs! 119831 : 
iStinebring et alJ[l990t iLommenl 120021 : 1 Jenet et al.ll2006f ) 
and the develop ment of improved analysis techniques for 
future searches dJenet et all 120051; lAnholm et all 120091: 
Ivan Haasteren et al.ll2009t) . A detection of the black hole 
stochastic background in the pulsar timing band may im- 
prove rate estimates for binary black hole mergers in the 
Laser Interferometer Space Antenna (LISA) frequency 
band (TWyithe fc Loebll2003T ). 

With any population there are always members that 
are nearer and brighter, and it has been predicted 
that several black hole binary systems should be indi- 



viduall y resolvable when the diffuse background is de- 
tected (|Sesana et al.ll2009l ). As a prelude to the first de- 
tection, upper bounds have been placed on the maximum 
amplitude of i ndividual systems us ing the Parkes Pulsar 
Timing Array (jYardlev et al. 2010;). When a detection is 
made, it has been shown that an analysis of the timing 
residuals imparted at the Earth can be used to constrain 
the amplitude of the signals to ~ 30%, and the direction 
to Afl ~ 40deg 2 for a 100 pulsar ar ray and a network 
signal to noise ratio of SNR = 10 (jSesana fc Vecchiol 
l2010al| bf) . But the timing residuals imparted at the Earth 
only tell part of the story. In addition to the disturbance 
at the Earth there is also the disturbance at the pulsar to 
consider. The pulsar component of the signal is usually 
ignored as it depends on the poorly constrained distance 
to each pulsar. In cross-correlation studies the pulsar 
terms average to zero as the projected distance to each 
pulsar is different, resulting in a different frequency and 
phase response to individual binaries. We show that it 
is possible to include the pulsar terms in the analysis by 
enlarging the parameter space to include the distance to 
each pulsar in the array. This doubles the signal power 
and allows the measurement of the mass and distance 
to the binary. As an added bonus, the pointing accu- 
racy improves by an order of magnitude, improving the 
prospects for finding electromagnetic counterparts. 

The outline of the paper is as follows: Section [2] de- 
scribes the response of the detector to a signal from a 
black hole binary system. An overview of the Bayesian 
inference techniques used to estimate the errors in the pa- 
rameter recovery is given in Section [3] Section [4] displays 
and discusses the results. The main results are summa- 
rized and future directions are outlined in Section [5] 

2. GRAVITATIONAL WAVES FROM 
SUPERMASSIVE BLACK HOLE BINARIES 

Assuming circular orbits, the orbital velocity of a black 
hole binary in the PTA frequency band scales as 
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where M = m\ + mi is the total mass, / is the grav- 
itational wave frequency, and we are using units where 
G = c = 1. We conclude from this that the orbital 
dynamics is only mildly relativistic, and that the grav- 
itational wave emission is well described by the lowest 
order post-Newtonian formulae. Another way of seeing 
this is to look at the time to merger, which scales as 
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Here M = (mima) 3 ' 5 /M 1 / 5 is the chirp mass, which 
ranges from 0.44M for equal mass systems to 0.22M for 
mass ratios of 1:10. These considerations suggest that 
the gravitational wave signal can be modeled as a simp le 
sinusoid of fixed fr equency (jSesana fc Vecchiol I2010al rbl: 
lYardlev et al.l [2010h . Allowing for moderate orbital ec- 
centricity introduces the complication of having to con- 
sider multiple sinusoids at harmonics of the orbital pe- 
riod, but the essential picture is unchanged. 

A more important effect that was not considered in 
these earlier studies is the contribution to the signal from 
gravitational waves disturbing the pulsars. This intro- 
duces a new time-scale into the problem in the form 
of the projected Earth-pulsar distance, d (1 — cos(fi)), 
where d is the distance to the pulsar and \x is the an- 
gle between the line of sight to the pulsar and the line 
of sight to the black hole binary. This time-scale is typ- 
ically of order a few thousand years. When the pulsar 
term is included, the effective baseline grows from tens, 
to tens of thousands of years (the temporal equivalent 
of aperture synthesis). The extended baseline makes it 
possible to measure the frequency change, and hence, the 
chirp m ass. The minimum dete ctable rate of frequency 
change (|Takahashi fc Setol l2002) 
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sets the minimum chirp mass that can be measured: 
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With the chirp mass determined, the amplitude of the 
signal can be used to solve for the distance to the black 
hole binary. The fi dependence in the effective distance 
provides additional information about the sky location of 
the binary, and this, combined with the increased SNR 
from using the full signal, leads to a significant improve- 
ment in the angular resolution of the PTA. 

The GW signal from a mildly relativistic black hole 
binary on a circular orbit is characterized by eight pa- 
rameters: the distance to the BH binary D; the sky lo- 
cation 4> and cos(0); the angular frequency lu = 7r/ of 
the binary orbit when observations begin at Earth; the 
orbital inclination cos(t); the orbital phase at the line of 
nodes 9 n ; the orientation of the line of nodes </>„ and the 
chirp mass A4 . The signal can be wr itten as the sum of 
two sub-signals ( Helli ngsl 1 1 98 ll [19821 ) . which we refer to 
as the pulsar signal and the Earth signal. The former is 
due to the disturbance caused by the gravitational wave 



at the pulsar, the later to the disturbance at the Earth: 

R(t e ) = r p (t p ) + r e (t e ). (5) 

In the plane wave limit, which pertains when D 3> d, t p 
is given by 

t p « t e - d (1 - cos(/x)) , (6) 

where t e is the time at the Earth and /i is the angle 
between the line of sight to the pulsar and the line of 
sight to the binar y. The two parts of the residuals ar e 
explicitly given in lWahlquistl (|1987f ): Uenet et all (|2004f > : 
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Here e + ' x are the GW polarization tensors, a is the unit 
vector pointing from the Earth to the pulsar, and 



with 



r+(t)=a(t) (A(t) cos(20) - B(t) sin(20)) 
r x (t) = a(t) (A(t) sm{2<t>) + B(t) cos(2<jf>)) 



A(t) = ~ sin [2 (Q(t) - 9 n )] [3 + cos(2t)] 



B(t) = 2 cos [2 (6(t) - 6„)] cos(t). 
The amplitude a(t) can be expressed as: 
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Finally, the orbital frequency and the orbital phase 
evolve according to 
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Fig. 1. — Residuals for a MBHB with chirp mass M = V2 X 
10 8 Mq, frequency of the gravitational perturbation at Earth f = 
5.0736 X 10 — 8 Hz and frequency of the gravitational perturbation 
at the pulsar / = 4.7145 X 10 -8 Hz 
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Unless the angle /i is very small, the signal at the pulsar 
will have a measurably different frequency compared to 
the signal at Earth, since the binary system will have 
evolved in the elapsed time At = d(l — cos(/i)) (Note: 
this quantity can be positive or negative depending on 
who "sees" the signal first). Figure [2] shows an example 
of the noise free timing residuals obtained from a MBHB 
system with A4 = v2 x 1O 8 M and initial gravitational 
wave frequency at the Earth f = 5.0736 x 10 _8 Hz. The 
system is located at a distance of D — 1.77 Mpc with 
a sky location of <f> = 0.3 and 6 = 1.42. The orbital 
inclination is l — tt/2, and the initial orbital phase and 
orientation of the line of nodes are both set at 7r/3. The 
superposition of the two signals is apparent from the beat 
envelope. It is possible to recover all the parameters 
characterizing the MBHB from the timing residuals seen 
in a PTA, so long as the pulsar contributions to the signal 
are used. 

At first sight it may seem impossible to include the pul- 
sar terms in a coherent analysis. To have phase errors in 
the pulsar terms less than AO radians requires that we 
know the distance to each pulsar to order Ac? ~ AG//. 
Setting A6 = 0.1 and / = 10~ 8 Hz gives Ad - 0.1 par- 
sees, or a fractional error of Ad/d ~ 0.01%. This is far 
smaller than the accuracy that is currently available from 
electromagnetic observations. Techniques such as paral- 
lax measurements and astrometry ach i eved a precision 
of about 10%. Recently I Verbiest et al.l ((200 8) have been 
able to estimate the distance to a specific pulsar with 
1% error using the kinematic distance derived from pul- 
sar timing data. This accuracy is not typical, and the 
method is only valid for nearby pulsars. But it turns out 
that highly accurate pulsar distance estimates are not 
required if we include the distance to each pulsar, di, as 
model parameters to be solved for from the GW data. 
The technique works as follows: for any one pulsar the 
phase matching of the pulsar terms produces a series of 
secondary maxima in di corresponding to 27r increments 
in the accumulated phase. As the estimate for di moves 
further away from the correct solution along this line of 
secondary maxima the predicted GW frequency and am- 
plitude of the pulsar term starts to deviate from the true 
value, and the height of the secondary maxima drop. In 
other words, it is the overall frequency /amplitude enve- 
lope that fixes the distance to the pulsar, and not the 
phase matching. Indeed, we will see that the secondary 
maxima are so close together that the individual peaks 
in the posterior distribution for di blend together into 
a continuum. The blending is even more pronounced in 
the marginalized posterior distributions for di where the 
correlations between di and /j. are integrated out. We 
will see that it is possible to use the GW data to provide 
estimates for the pulsar distances, while simultaneously 
deriving useful estimates for the chirp mass and black 
hole location. 

It was recently pointed out by iDeng fc Finn! (|2010l ) 
that the timing residuals may show departures from the 
plane wave approximation used to derive equation ([8]). 
They showed that the curvature of the wavefronts in- 
troduces parallax effects which can be used to measure 
the distance to the black hole binary and the pulsars in 
the array, independent of including the pulsar term in 
the signal. Accounting for the wavefront curvature will 



improve our estimates, which neglect this effect. 

The common challenge with all GW detectors is to 
isolate a signal which is relatively weak compared to 
the noise surrounding it. The noise is composed of in- 
strumentation or measurement noise and the background 
confusion noise generated b y the superposition of all the 



unresolved MBHBs signals (iRaiagopal fc Romanil 
Phinnevl I200H; Uaffe fc Backed 120031: Uenet et all 



1995 



2005 



2006: ISesana et all 120081) 7" The measurement noise is 



simulated as a white noise with a standard deviation 
of a = 100 ns for all pulsars. An array of 20 pulsars 
each with a timing rms of about 100 ns is a reasonable 
estimat e of what will be achieved by PTAs i n the near 
future (|Manchesterll2008l: iVerbiest et al.ll2009f ). It is also 
thought to be close to t he threshold for th e detection of 
the diffuse backg round dJenet et all 120051). Models for 
the ba ckground (jSesana et al.1 120081 : ISesana fc Vecchid 
l2010al ) suggest that its power spectral density can be 
described as: 
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where /i*, /* and 7 are model dependen t . Th is result 
differs from earlier studies (e.g. P hinnevl (|200lD ) which 
predicted a simple f~~s power-law. In Figure [2] we 
compare the power spectral density of the signal gen- 
erated by the MBHB described above with white instru- 
ment noise and four different models of background noise 
(e.g. VHM, VHMho pk, KBD, BVRhf for a description 
of th e se models see (ISesana et al.l l2007t IVolonteri et all 
20031: iKoushiappas et all 12004 iBegelman et all 120061: 
Volonteri et alll2006D l 
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Fig. 2. — Power spectral density of the full (solid line) and trun- 
cated (dashed line) signal from a binary compared with the power 
spectral density of the instrument white noise and four models of 
background noise (VHM, VHMhopk, KBD, BVRhf). The signal 
power is shown for a single pulsars in the array. The cumulative 
SNR (over the full array) is 20 for the full signal and 14 for the 
truncated signal 

The square of the signal to noise ratio doubles when 
using the full signal compared to the truncated signal, 
which signifies that including the distances to the pulsars 
in the parameter space could contribute to the detection 



4 



of more sources. For most of the frequency window of 
interest, the background no ise is smaller than both the 
signal and the white n oise (jRaiagopal fc Romanil 119951 : 
Uenet et all I2005L [2006). For the sake of simplicity, we 
will ignore the background for now, as it will not greatly 
affect the estimation of the parameters for a detected 
source. One would need to include it however in order 
to determine the number of detectable sources. 

3. BAYESIAN INFERENCE AND PARAMETER 
ESTIMATION 

To determine how well we can estimate the parame- 
ters of a MBHB from noisy pulsar timing data we need to 
compute the posterior distribution for the model parame- 
ters. The posterior distribution p(x\s) gives the probabil- 
ity of observing parameters x given data s. It is defined 

by 

p{s\x)p{x) 
J p{s\x)p(x)dx 

Here p{x) is the prior distribution, which is the mathe- 
matical representation of the prior knowledge of the sys- 
tem, and p(s\x) is the likelihood evaluated at x. It is the 
probability that a system with parameter x will yield a 
signal s in the detector. 

To generate the posterior distribution , a Markov Chain 
Mont e Carlo ( M CMC ) algorithm (Metropolis et all 
1951 lHastingsl [l970i: iNewman fc Barkemal I1999T 
van Haasteren et al.ll2009D is used to explore the param- 
eter space. The MCMC consists of proposing "jumps" 
from one location Xi in the parameter space to another 
Xi + ±. The jumps have a finite probability K,(Xi-\-i\xi) of 
being accepted that is given by the Hasting ratio: 



of Gaussians centered around the measured value df M 
with a standard deviation <x; = adf M : 



where 



H = 



K(x i+ i\xi) = min(l, H) 

p{x l+ i)p(s\x i+1 )q(xi\x i+1 ) 
p(xi)p(s\xi)q(x i+ i\xi) 



(17) 
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Here q{y\x) is the proposal distribution: the probability 
that a jump from x to y will be proposed. For a more 
detail ed description of MCMC techniques see lGamermanl 
(1997ft . 

3.1. Prior distribution 

Some information is known prior to the analysis. For 
instance, there is a higher probability the source will be 
located far from the Earth, as the area of a sphere in- 
creases as the square of its radius. To account for this, 
the prior distribution on D is chosen to be proportional 
to D 2 out to some large distance -D max that is well outside 
the horizon for PTA detections. The sources are assumed 
to be uniformly distributed on the sky, with random ori- 
entations. A more informative prior on the distance and 
sky location could be built using galaxy catalogs. For the 
distances to the N pulsars, the prior follows from electro- 
magnetic observation^]], which we take to be a collection 

1 Paul Demorest has pointed out to us that for distance estimates 
derived from parallax measurements, it would be better to choose 
a prior that is Gaussian in the parallax 1/di. However, if the 
fractional error in di or 1/di is small, the difference between using 
distance and parallax to define the prior will also be small. For 
a recent discussion o f the statistics of parallax measurements see 
IVerbiest efaT] 112010) . 
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where di is the proposed distance to the pulsar i. For our 
simulated PTA we draw d from the range 0.5 — 1.5 kpc, 
and include errors in the estimated d EM that are con- 
sistent with our choice of prior. We will consider some 
examples where the fractional distanced error a — 0.1 for 
all the pulsars, and other examples where a takes differ- 
ent values for different pulsars. We will focus on arrays 
with N = 20 pulsars. 

3.2. Likelihood 

The likelihood p(s\x) is by definition the probability 
of observing a signal s from a source with parameters x. 
For Gaussian noise it is given by: 



p{s\x) = C exp 



±((s-R(x))\(s-R(x))) 



, (20) 



where C is a normalization constant, R(x) is the wave- 
form described in Section[2]for a set of parameters x, and 
the noise weighted inner product is defined as: 

(a\b) = 2 [ X (a* (/)&,(/) + <k(f)b*Ajj) S m H.f)df. 



(21) 

A summation over the indexes i and j from 1 to N is 
implied. S ni j(f) is the spectral density of the noise cor- 
relation matrix CVj(r) — J ni(t)rij(t + r)dt where rii(t) 
is the noise in the signal from pulsar i. S ni j is not typ- 
ically diagonal since the background noise is correlated 
between the pulsars. However the instrument noise is un- 
correlated. Since we ignore the background noise, S ni j 
becomes diagonal. 

3.3. Proposal distribution 

A combination of six proposal distributions q{y\x) is 
used. The first makes use of the Fisher Information Ma- 
trix approximation to the posterior distribution. The 
Fisher matrix indicates the level of correlation between 
the parameters, and the diagonal elements of its inverse 
give a rough approximation of the error expected in the 
estimation of each parameter. The Fisher matrix is the 
expectation value of the negative Hessian of the log pos- 
terior evaluated at the posterior mode: 

r,j = -(dtdj log p (x\s)) = (dtRldjR) - didj log p (x) 

In the current setting the prior only contributes to the 
diagonal elements of the Fisher matrix, adding a term 
1/af to the elements representing the pulsars distances. 
The eigenvectors of the Fisher matrix define a new set 
of uncorrelated parameters. The eigenvalues indicate the 
curvature of the likelihood surface along each eigenvec- 
tor. When the curvature is high, the likelihood changes 
a lot for a small variation of the "eigen-parameter" , and 
a big jump is unlikely to get accepted. To effectively 
explore the likelihood surface, jumps along those eigen- 
vectors are proposed, which are scaled by their eigen- 
values. The second proposal distribution is similar, but 
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with a bigger scaling. Since the jumps are bigger they 
are less likely to be accepted, but when accepted they 
explore the parameter space faster. The third proposal 
consists of drawing a new parameter set from the prior 
distributions. The fourth proposal pick selects one of the 
pulsars and draws a new pulsar distance from the prior 
distribution. The fifth proposal distribution is similar 
to the fourth, but uses a proposal centered around the 
true value. It prevents the MCMC from searching exclu- 
sively around the values predicted by the prior. Finally, 
tiny jumps are sometimes proposed along each parame- 
ter. These tiny jumps are very likely to be accepted, and 
this helps prevent the chain from getting stuck in a loca- 
tion where the Fisher matrix estimates are particularly 
poor approximations to the posterior distribution. 

3.4. Parallel tempering 

A parallel tempering scheme (jSwendsen fc Waiig1ll986( ) 
is implemented to improve mixing and convergence of the 
Markov Chains. It consists of running a number of chains 
in parallel, each of them with a "temperature" T = 1//3 
which modifies the likelihood: 

p^{s\x^)=p{s\xf. (23) 

The temperature effectively "smoothes" the likelihood 
map. The higher the temperature, the more likely a jump 
is going to be accepted. The chains then communicate by 
swapping with each other. The swaps have a probability 
of being accepted given by a Hasting ratio: 

H = p(s|xi,/3j)p(s|xj,A) 
p(s\&i,Pi)p(s\xj,f)j) ' 

The indices i and j refer to the individual chains. Here 10 
chains are used, with temperatures exponentially spaced 
between T = 1 and T = 10. 

4. RESULTS 

The addition of the pulsar distances to the parameter 
space allows us to use the pulsar signals in the analy- 
sis. It follows that information about the evolution of 
the orbital frequency over a long period of time and the 
relative phases between the Earth signal and the pulsar 
signals can be extracted. This information is enough to 
independently determine the distance to the binary and 
the chirp mass. 

We use a simulated PTA comprised of N — 20 ran- 
domly chosen pulsars drawn uniformly in sky location 
and with distances in the range 0.5 — 1.5 kpc. The pa- 
rameters describing the array are listed in Table [TJ With 
20 uniformly distributed pulsars the array has a fairly 
uniform antenna pattern. Increasing the size of the ar- 
ray to 40 pulsars would further improve the sky coverage, 
but the effect is not very significant. To investigate these 
effects we studied the the spread in SNR across the sky 
for a particular source, by considering 3072 different sky 
locations using 50 randomly drawn arrays of pulsars. In 
each case the SNR corresponding to the sky location was 
calculated, and normalized by the average SNR to pro- 
duce the histograms seen in Figure [3l The variation in 
SNR across the sky is smaller for the larger PTA, but not 
by a significant amount. In either case, we do no expect 
to see a significant correlation between the sky location 
and black hole distance parameters, and this expectation 



is borne out in our analysis. The distance and sky loca- 
tion of the 20 pulsars used in our analysis are listed in 
Table [U 
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FlG. 3. — Typical SNR spread across the sky for MBHB signals 
observed by arrays of 20 and 40 pulsars 





d (kpc) 


6,em (kpc) 




cos 


pulsar 1 


0.987 


1.010 


1.17 


-0.27 


pulsar 2 


1.391 


1.371 


5.49 


0.95 


pulsar 3 


1.008 


1.299 


6.28 


-0.25 


pulsar 4 


1.128 


1.057 


5.21 


-0.99 


pulsar 5 


1.075 


1.111 


5.87 


0.63 


pulsar 6 


0.651 


0.629 


5.23 


-0.48 


pulsar 7 


0.957 


0.895 


0.02 


0.67 


pulsar 8 


1.138 


0.955 


5.13 


-0.37 


pulsar 9 


1.133 


0.995 


3.24 


0.10 


pulsar 10 


0.516 


0.771 


5.23 


-0.72 


pulsar 11 


1.184 


1.084 


4.53 


0.09 


pulsar 12 


1.674 


1.491 


5.85 


-0.52 


pulsar 13 


1.455 


1.437 


3.57 


0.48 


pulsar 14 


0.544 


0.507 


2.83 


-0.70 


pulsar 15 


0.897 


0.925 


0.34 


-0.28 


pulsar 16 


0.756 


0.688 


4.71 


-0.98 


pulsar 17 


1.318 


1.331 


6.21 


0.09 


pulsar 18 


0.882 


0.920 


1.46 


0.14 


pulsar 19 


0.676 


0.641 


0.23 


-0.64 


pulsar 20 


0.852 


0.757 


5.60 


0.09 



TABLE 1 

The position and distances for the 20 pulsars that define 
the array used in the analysis. here d is the true distance 
to the pulsar, while &em is the prior estimate of the 
distance from traditional astronomical methods. 



The distance to the binary is only present in the am- 
plitude of the signal. It makes its determination very 
dependent on the inclination angle i. For values of the 
inclination angle close to or tt, the determination of the 
distance to the black hole binary system is very poor. At 
t = and t = 7T, the Fisher matrix becomes singular. 
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Figure HI displays the error in the distance to the binary 
predicted by the Fisher information matrix as a function 
of the inclination. 
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Fig. 4. — Error in the determination of the distance to the binary 
as a function of the inclination. The error is estimated using the 
Fisher information matrix. The breaking of reflection symmetry 
about l = 90° comes from using a fixed angle to the line of nodes, 
which means that polarization angle changes with t. 

For an ideal value of 7r/2 for the inclination and a 
SNR of 20, our MCMC predicted an error of - 7% 
in the distance estimation, which is in agreement with 
the Fisher prediction. The error on the chirp mass was 
~ 2%. The orbital frequency was extremely well de- 
termined (~ 0.3%). Figure [|] shows the marginalized 
posterior distribution for the black hole distance and the 
chirp mass for three binaries with different chirp masses, 
1O 8 M ,5 x 1O 8 M and 1O 9 M respectively (Sources 1, 
2 and 3 in Table . The marginalized posterior distri- 
butions are represented by the boxed histograms while 
the smooth line represents the Fisher matrix estimates. 
The distances were chosen in order to ensure a SNR of 
20 for each case. The heaviest source is sufficiently rela- 
tivistic that higher order post-Newtonian effects may be 
detectable, but we defer this analysis to a future study. 

The Fisher information matrix is seen to provide a 
good approximation to the posterior distribution for the 
black hole parameters. The differences can in part be 
attributed to imperfect convergence in the MCMC runs, 
which suffer from the high dimensionality of the param- 
eter space and strong correlations between many of the 
model parameters. In some cases the posterior distribu- 
tion peaks are shifted from their injected values. This 
is caused by the priors for the distances to the pulsars, 
which are not always peaked close to their true values, 
as we will show below. 

Unfortunately, most black hole binaries will not have 
an inclination of i = tt/2. For a more realistic error pre- 
diction, we perform the same analysis for similar binaries 
whose inclinations are this time chosen to be 7r/4. We 
otherwise use the same sets of parameters (Sources 4, 5 
and 6 in Table [2]) . The distances D change slightly in 
order to conserve an SNR of 20. The results are given in 
Figure 0] 

The error for the distance D is predictably bigger. It 
rises from less than 10% to about 15% The Fisher ma- 
trix estimates the error to be as large as 30%. Again one 
has to take into consideration the Fisher matrix gives 
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Fig. 5. — Posterior distributions of D and Ai with i = 7r/2 for 
three binaries with different chirp mass: 10 s Mq (first row), 5 X 
10 8 M Q (second row), 10 9 M Q (third row). The distance to the 
three binary are normalized so that the three SNR are equal to 20. 
The MCMC derived posterior distributions (boxes) are compared 
to the Fisher matrix estimates (dashed line). 



only a rough approximation to the posterior distribution. 
The posterior distributions obtained for the three differ- 
ent binaries are consistent, which indicates the MCMCs 
have reasonably converged. The measurements of the 
frequency and chirp mass were not significantly affected 
by the change in inclination, which confirms they are 
weakly correlated. 





D (10* kpc) 


M (10" Mq) 


woO-Q-'s- 1 ) 












1 


1.05 


1.0 


1.328 












2 


16.ti 


3.0 


1.328 












3 




10 


1.328 


tt/2 






1.05 


6.08 


4 


2.37 


1.0 


1.328 


tt/4 






1.05 


6.08 




36.2 


5.0 


1.328 


*/4 


3.98 


U.*:,070 


1.05 


6.08 


C 


121 


10 


1.328 


7T/4 


3.98 


-0.83070 


1.05 


6.08 



TABLE 2 

The parameters for the 6 black hole binary examples we 

are considering. the merger times t c can be computed 
from these parameters: sources 1,4 tc = 4.39 x 10 4 years; 
Sources 2,5 t c = 3.0 x 10 3 years; Sources 3,6 t c = 9.46 x 10 s 

YEARS. 



In addition to the decoupling of the binary distance 
and chirp mass from the amplitude, the addition of the 
pulsar terms increase considerably the precision of the 
determination of the binary position. The pulsar terms 
are evaluated at t p , which is the time at which the distur- 
bance from the gravitational wave occurred. It is given 
by t p = t — d(l — cos(/i)), where (i is the angle between 
the line of sight to the pulsar and the line of sight to the 




D(lQ 3 kpcJ MflQ a Wg) 



Fig. 6. — Posterior distributions of D and Ai with t = 7r/4 for 
three binaries with different chirp mass: 10 s Mq (first row), 5 X 
10 s Mq (second row), 10 9 M Q (third row). The distance to the 
three binary are normalized so that the three SNR are equal to 20. 
The MCMC derived posterior distributions (boxes) are compared 
to the Fisher matrix estimates (dashed line). 



binary. It is therefore a function of the position of the 
binary. The pulsar signals give previously non-existent 
information on the chirp mass and distances to the pul- 
sars, but also help refine the measurement of the binary 
sky location. To illustrate this, Figure [7] compare the 
posterior distribution with the Fisher estimation for the 
two position parameters, for a case in which the full sig- 
nal was used and a case in which the signal from the 
disturbance at the pulsar was omitted. When the pul- 
sar signals are added, the SNR naturally increases. Here 
the distance to the binary was normalized each time to 
get a SNR of 20 for both cases. The effect from the 
information encoded in the new term in the signal can 
then be dissociated from the effect due to the increase 
in the SNR. A mid-range chirp mass and frequency were 
used (M — 5 x 1O 8 M and uj = 2ir/ 1.5 years), and an 
inclination of l = 7r/4. 

To compare our results to iSesana fc Vecchio (2010b), 
we calculate the pulsar timing array angular resolution 
given by 

Ail = 2vry / (Acos6'A0) 2 - (C"*> cose ) 2 , (25) 

where cosd is a off-diagonal term of the covariance 
matrix (the inverse of the Fisher information matrix). 
For the errors in the individual sky location parameters 
we use the Fisher matrix estimates, which are in very 
good agreement with the posterior distributions. For the 



Fig. 7. — Position errors for the full signal (upper panel) and the 
truncated signal (lower panel). 



truncated signal we find ACl = 41.3 deg 2 , which is consis- 
tent with[Sesana fc Vecchio! (|2010b| ). while the full signal 
(with the distance changed to preserve the same SNR) 
yields Afl — 6.5 deg 2 . The measurement is therefore im- 
proved by a factor of 6.4. Keeping the distance fixed and 
accounting for the increase in the SNR when the pul- 
sar terms is included, the angular resolution goes from 
An = 41.3 deg 2 to An = 2.6 deg 2 . This is a consider- 
able improvement (over an order of magnitude), which 
highlights the importance of utilizing the full signal. 

It was previously mentioned that the full signal also 
gives information about the distance to the pulsars. In 
the simulations described above we assumed a 10% error 
in each of the pulsar distances, which correspond to the 
high end of today's accuracy in measurements using par- 
allax and other methods. We expect the pulsar distances 
to be further constrained by the gravitational wave anal- 
ysis. The Fisher information matrix predicts that for 
most pulsars, the error in the distance measurement will 
be constrained to a few percent. The level of improve- 
ment depends on the location of the pulsar with respect 
to the line of sight to the binary. If the Earth , pulsar 
and binary are aligned, then t p — t, and no new informa- 
tion is given by the pulsar term in the signal. Figure [8] 
compares the Fisher matrix estimates to the MCMC de- 
rived posterior distribution and to the prior distribution, 
for a few relevant pulsars: 

A few things are clear at first sight. First, when the 
pulsar's sky location is close to the line of sight to the 
binary (cos fi « ±1), the prior is recovered from the pos- 
terior distribution, which means that no new information 
was acquired from the gravitational wave signal, as ex- 
pected. Then, for the pulsars that are not close to the 
line of sight, the peak of the posterior distribution can be 
shifted from the true value of the distance. This occurs 
when the true value (value at which the Gaussian ex- 
tracted from the Fisher information matrix is centered) 
is located a few sigmas away from the prior-predicted 
value. As mentioned earlier, these shifts will induce a 
shift in the posterior distribution of the binary distance, 
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Fig. 8. — Marginalized posterior distributions for the distances 
to 8 of the 20 pulsars in the array (boxes) compared to the prior 
distributions (dotted line) and the Fisher matrix estimates (dashed 
line). The distance priors have a standard deviation of 10%, which 
correspond to the confidence in the measurement from electromag- 
netic astronomy. For some of the pulsars, the gravitational wave 
signal slightly improves the distance determination. 



chirp mass and orbital frequency with respect to their 
true values. As a sanity test, the priors were centered to 
the right values for all the pulsars' distance. As antic- 
ipated, the posterior distribution for the binary param- 
eters are in this case centered on their right values as 
well. Finally, and maybe most importantly, the poste- 
rior distribution from the MCMC is consistently much 
broader than the error predicted by the Fisher informa- 
tion matrix. For most pulsars, the posterior distribution 
is as wide as the prior, which means that only limited 
information concerning the distances could be extracted 
from the gravitational wave signal. To explain the dis- 
crepancy, the pulsar term of the residuals is rewritten 
as: 

Mi 

r P^ = n cos M*p)*p + $) , (26) 

Duj a (tp) 

with t p — t — d(l — cos (j,). If one varies the distance to 



the pulsar d such that 

u(t'p)t' p = Lo(tp)t p + 2tt, (27) 

then only the amplitude and frequency are slightly 
changed. This results in a series of secondary maxima 
in the likelihood, which for this particular pulsar are 
spaced by Ad ~ 0.017 kpc, which correspond to 1.7% 
of the distance. The change in frequency and amplitude 
between adjacent maxima is small compared to the mea- 
surement uncertainty in these quantities, and it is not un- 
til multiple secondary maxima have been traversed that 
the likelihood drops significantly. The separation of the 
secondary maxima is comparable to the Fisher matrix 
prediction for the error in d, which means that where 
the Fisher matrix predicts the posterior distribution to 
drop, the MCMC will find another maximum, almost as 
good as the primary. The error in the detection of each 
local maximum blends with the error in the detection 
of its neighbor. They form an "error envelop" which is 
limited in its width by the change in the frequency and 
amplitude of the pulsar signal. Even though the gain in 
the precision of the measurement of the distances to the 
pulsars is not as significant as predicted by the Fisher 
matrix analysis, it is still noticeable for a few pulsars. 

To explore the limitation of the determination of the 
distances due to the periodicity of the residuals, we set 
the priors for five pulsars to be Gaussians with stan- 
dard deviations randomly drawn from the range a — 
[0.005,0.03]. For ten pulsars, the standard deviation is 
drawn from the interval a = [0.009,0.15], and for the 
remaining five pulsars, the prior is chosen to be less con- 
straining, a = [0.20,0.25]. Table [3] lists the estimated 
distances to the pulsars d l EM , their true distances d l , and 
the estimated error in the existing measurement repre- 
sented by the standard deviation of the prior distribution 

^prior ■ 



Pulsar 


d (kpc) 


Aem (kpc) 


ffprior 


1 


1.009 


1.010 


0.53 


2 


1.374 


1.371 


1.22 


3 


1.242 


1.299 


1.96 


4 


1.070 


1.057 


1.75 


5 


1.109 


1.111 


0.52 


6 


0.655 


0.629 


11.7 


7 


0.965 


0.895 


11.2 


8 


1.139 


0.955 


10.0 


9 


1.135 


0.995 


10.1 


10 


0.494 


0.771 


10.9 


11 


1.184 


1.084 


10.0 


12 


1.700 


1.491 


11.4 


13 


1.455 


1.437 


10.0 


14 


0.547 


0.507 


10.8 


15 


0.895 


0.925 


10.6 


16 


0.827 


0.688 


20.4 


17 


1.301 


1.331 


23.3 


18 


0.841 


0.920 


20.6 


19 


0.725 


0.641 


24.1 


20 


0.935 


0.757 


18.8 



TABLE 3 

The distances to the 20 pulsars in the array. Here d is the 

TRUE DISTANCE TO THE PULSAR, WHILE d EM IS THE PRIOR 
ESTIMATE OF THE DISTANCE BY TRADITIONAL ASTRONOMICAL 
METHODS, AND <J pT i OT IS THE ESTIMATED ERROR IN THE 
MEASUREMENT OF THE DISTANCE (AS A PERCENTAGE OF d EM ). 



The estimation of the pulsar distance is improved sig- 
nificantly for pulsars 7, 17, 18, 19, 20. The improvement 
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is largest for pulsars that were poorly constrained by elec- 
tromagnetic measurements, and did not lie close to the 
line of sight to the black hole (cos /i ~ 1). When the prior 
distance is far from the real value, the gravitational wave 
data picks up on the discrepancy, and the posterior distri- 
bution becomes double peaked (pulsar 10). The posterior 
distribution of pulsar 7 is identical whether the five first 
pulsars have a tight prior or not. This seems to indicate 
that the improvement on the determination of the dis- 
tances is still limited by the periodicity of the strength of 
the residuals with respect to the pulsars distances. For 
this reason, the posterior distribution remains broader 
than the Fisher estimation. Figure |H] display the poste- 
rior distribution for the distances to a few relevant pul- 
sars against the Fisher estimations and the priors. 




Pulsar IS (kpe), cos(u) = -0.57 Pulsar 19 |kpc), cos(u) = 0.18 



Fig. 9. — Marginalized posterior distributions for the distances to 
8 of the 20 pulsars in the array (boxes) compared to the prior dis- 
tributions (dotted lines) and the Fisher matrix estimates (dashed 
line). Pulsars 1 — > 5 were assumed to have distances that were well 
determined by electromagnetic observations. As a consequence, the 
measurement of the distances to the some of the remaining pulsars 
can be significantly improved by folding in the gravitational wave 
analysis of the timing residuals. 



The errors on the binary parameters 
(D, A4, w, (j>, cos 9) were not significantly affected 
by the change in the pulsars priors. 



5. CONCLUSIONS 

We have presented a novel method of analyzing binary 
black hole signals using pulsar timing data. By including 
the distances to the pulsars as model parameters we are 
able to incorporate the "pulsar term" in the gravitational 
wave signal, which allows us to detect the decay of the 
orbit, and hence the chirp mass. This in turn allows us to 
convert the amplitude of the signal into a measurement 
of the distance to the black hole binary. Including the 
pulsar terms doubles the signal power and reduces the 
pointing error by an order of magnitude. For detections 
with a network SNR = 20, it should be possible to mea- 
sure the distance to < 20%, the chirp mass to < 5%, and 
the sky location to Af2 < 3deg 2 . In something of a role 
reversal, the gravitational wave observations can improve 
the distance estimates to pulsars in the array. The im- 
provement can be significant for pulsars whose distances 
are originally poorly estimated. It also follows that the 
more MBHBs detected, the more trustworthy the estima- 
tions for the distances to the pulsars are, which in turns 
allow for stronger constraints on the MBHBs detected. 

Topics for future research include the study of sys- 
tems with eccentric orbits, and the possibility of detect- 
ing higher order post-Newtonian effects (which encode 
information about the mass ratio and spin of the black 
holes) for higher frequency or more massive systems. It 
would be interedting t o see how includin g the effects 
of wavefront curvature (jDeng fc Finn 2010) improve the 
distance estimates. Future studies could also consider the 
effect of the confusion noise from the unresolved MBHB 
background, which will ultimately determine the number 
of individual systems that can be resolved. 
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